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A spin-orbit-coupled Bose-Einstein-condensed cloud of atoms confined in an annular trapping 
potential shows a variety of phases that we investigate in the present study. Starting with the non¬ 
interacting problem, the homogeneous phase that is present in an untrapped system is replaced by a 
sinusoidal density variation in the limit of a very narrow annulus. In the case of an untrapped system 
there is another phase with a striped-like density distribution, and its counterpart is also found in 
the limit of a very narrow annulus. As the width of the annulus increases, this picture persists 
qualitatively. Depending on the relative strength between the inter- and the intra-components, 
interactions either favor the striped phase, or suppress it, in which case either a homogeneous, or a 
sinusoidal-like phase appears. Interactions also give rise to novel solutions with a nonzero circulation. 

PACS numbers: 05.30.Jp, 03.75.Mn, 67.85.Fg, 67.85.Jk 


I. INTRODUCTION 


The effects associated with spin-orbit coupling have 
long been known and studied in various physical systems, 
includiM atoms, solids, quantum dots, etc. [lH3! ^^.d 
nuclei |8|. Spin-orbit coupling often plays an important 
role in semiconductor nanostructures, where the field of 
spintronics has led to many new applications Q. Con¬ 
cepts known from semiconductor nanoscience are nowa¬ 
days often found to be transferable to confined atomic 
quantum gases, with an ever increasing interest in the 
properties of “atomtronic” devices that make use of cold 
atoms in a similar way like electrons in solids [iMi. 
Semiconductor “spintronic” devices are to a large extent 
controlled by the Rashba or Dresselhaus spin-orbit cou¬ 
pling in the heterostructure, giving rise to a number of 
intriguing spin-dependent transport phenomena M- By 
means of laser-coupling techniques an analogue to spin¬ 
tronics has now become possible with the achievement of 
artificially-induced spin-orbit coupling in ultracold quan¬ 
tum gases of neutral atoms with the exciting possi¬ 
bility of creating “atom-spintronic” devices that may also 
make use of the superffuid properties of quantum gases. 
Recently, such systems have been experimentally real¬ 
ized in Bose as well as in Fermi 231 gases. Their 

unique properties and their fundamental and application- 
related prospects have inspired intense theoretical efforts 
[234^ . Furthermore, in a series of other experiments, 
it has become possible to trap atoms in annular/toroidal 
traps, see, e.g., Refs, 


Motivated by the above experiments, we investigate 
the lowest-energy states that appear in a spin-orbit- 
coupled Bose-Einstein-condensed cloud of atoms that is 
trapped in an annular potential, using the mean-field ap¬ 
proximation. A lot of work has been done in homoge¬ 
neous [ 13 , 113 , IsB - lssjl and in harmonically-trapped sys¬ 


tems m, m, m, m, no, 113, m-im • in the presence of 
an annular potential, however, this problem becomes sur¬ 
prisingly challenging. As we show below, even in the ab¬ 
sence of interactions the eigenvalue problem due to the 
spin-orbit coupling is non-trivial. 

We proceed in three steps: First of all, we ignore the 
interactions and consider the case of a very tight annulus, 
developing an effectively one-dimensional theory. Within 
this model we integrate over the transverse degrees of 
freedom, which not only simplifies the problem numer¬ 
ically, but it also gives some additional insight into the 
properties of the system. Furthermore (Sec. Ill), the re¬ 
sulting equations have some limiting analytic solutions. 
We then solve the eigenvalue problem in an annulus of a 
finite width (Sec. IV). The interactions are incorporated 
in the last step (Sec. V). We pay special attention to the 
interplay between the (single-particle) effect of spin-orbit 
coupling and the atom-atom interactions; we stress that 
the parameters that are associated with both of them 
are realistic and controllable experimentally in atomic 
systems. 


II. MODEL AND METHOD 

We consider a Bose-Einstein-condensed cloud of atoms, 
confined to two dimensions. Thus, the single-particle 
Hamiltonian of the system is 

N 

[-?l"V?/(2M) + V{ri) + Vso(d)] • (1) 

i=l 

The external confining potential has the form of an an¬ 
nular potential, V{r) = Muj^{p — R)'^j 2 , where p is the 
radial coordinate, M is the atom mass, oj is the frequency 
of the potential, and R is the mean radius of the annu- 
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lus. The term Vso associated with the spin-orbit cou¬ 
pling that we consider here is the one that corresponds 
to the experiments of the Spielman group [^ . In this 
case the so-called Rashba-like and Dresselhaus-like spin- 
orbit contributions have an equal strength, and thus Vso 
takes the one-dimensional form 


;;2jL2 

Vso{t^ = ^ 
sov y 2M 


+ +-h6ay. 


( 2 ) 


the population imbalance between the two components. 
The most “interesting” is the spatially-dependent term, 
which is proportional to fco ■ It is this term that gives rise 
to the striped phase. 


III. NON-INTERACTING PROBLEM IN THE 
LIMIT OF QUASI-ONE-DIMENSIONAL MOTION 


Here fcg is the wavenumber of the Raman laser beams, 
Px is the linear momentum of the atoms in the x direc¬ 
tion, fj-l (with j = X, y, z) are the Pauli matrices, O is the 
Rabi frequency and S is the detuning. The experimen¬ 
tal studies of spin-orbit coupling generally choose small 
values of 6 [U, and thus we ignore the corresponding 
term in what follows below. Furthermore, we here orient 
the spin axes as in Ref. 22|, but most recent theoreti¬ 
cal work uses a cyclic spin rotation ay ^ az, etc. 

Therefore, the system is described as an effective pseudo¬ 
spin- 1 / 2 , and the condensate wave function is written as 
a two-component vector [$^,< 1 )^]^ (“up” and “down”). 

The interactions are modelled via a short-ranged s- 
wave potential of the form 


Mn. = f E 


dxdy + da;d?/. (3) 


While in general there are three different interaction 
strengths between the “up” and “down” components, 
t^ 2 e present problem we set = 
9 ii = 9 tune the ratio 9 / 9^1 (both assumed to be 
positive). Our total Hamiltonian is thus H — Hq + Vint, 
which yields two coupled Gross-Pitaevskii-like equations 
for the two components of the order parameter. 


' Hn 

H 12 ' 


■ $4, ■ 

= P 

■ $4, ■ 

TJ * 

. -^12 

H22_ 





where Hu = It?/{2M) -|- p^/(2M) + V{p) + Ml/2 + 
H 22 = T?kl/{2M)+p^/{2M) + V{p)- 
hn/2+g\^??+g^l\^^?, Hi 2 = -ihkoPx/M-ih5/2, and 
p is the chemical potential, which is determined from the 
condition that the total occupancy of the two pseudospin 
components and iVj, is equal to the total number of 
atoms N. 

For the specific form of the spin-orbit coupling consid¬ 
ered and in the absence of any trapping potential and of 
interparticle interactions, the density distribution of the 

S 5 may be either striped, or homogeneous [Ullllli- 
. This is determined by the ratio between hil and 
the recoil energy En = t?kQ/{2M) If the ratio 

Ml/{AE}i) is smaller or larger than unity, the system 
is in the striped, or in the homogeneous phase, respec¬ 
tively. Furthermore, the expectation value of Vso in 
some state [$- 1 -, $ 4 .]^ is Eso — E^kQ/{ 2 M) = h^/ 2 {N-f — 

N^)-{h‘^ko/M) J (^<^*d<i>:^/dx - ^ld<^^/dx^ dxdy. The 

term proportional to kg is constant. The term propor¬ 
tional to n is spatially-independent and is responsible for 


For very strong transverse confinement, the atoms re¬ 
side in the ground state associated with the transverse 
degrees of freedom for motion along the annulus. This 
fact allows us to make an ansatz for $ 4 -, 

(5) 

with a Gaussian density distribution in the trans¬ 
verse direction, i.e., no{p) = \4>o(p)? — exp[—(p — 
i?)^/ag]/(i/^aoi?), where ap is the oscillator length, oq = 
(h/MurY''^. Integrating over the transverse degrees of 
freedom and assuming that i? ^ oq we thus develop the 
following quasi-one-dimensional eigenvalue problem for 
^' 4 - and 'I'j, (i.e., two coupled differential equations). 


de^ ^ 2^^ + 


h^ko 

MR 


1 „ ■ n 9 

-cos9 + sm9— 

2 o9 






1 , 


2MR'^d9^ 2^^ 


t? ko 

MR 


1 d 
-cos9 + sm9— 

2 o9 


= 0 , 


( 6 ) 


where E is the eigenenergy. 

In the absence of interactions there are three en¬ 
ergy scales, f?/{MRY, f?ko/{MR), and h^. The 
above coupled equations have an analytic solution in 
the limit where ??ko/{MR) is much smaller either than 
T?/{MR?), or than fiVt (we stress that for fcp = 0 the two 
equations decouple, while for “small” values of kgR they 
couple perturbatively). Defining the two dimensionless 
parameters koR ^ I (i.e., the ratio between the second 
scale and the first), and A = hko/(MilR) (i.e., the ratio 
between the second scale and the third), we find that for 
koR < I, 


VI/ 4 , 


koR 1 
y/^l + e 


cos 9, 



1 — -(fcoi?)^(l -I- e) cos20 


(7) 


The eigenenergy is E/{hfl) = —1/2 — (/coi?)^/[2e(l -I- e)], 
where e = 2MR^Vl/Ti is defined as the ratio between Ml 
and T?/{2MRY- If a = hfl/En, then e = a(fcoi?)^. For 
typical values of a of order unity, e ~ (koR)'^. A second 
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analytic solution is found for A <C 1, 


vft = - 


= 


A 




cos 9, 




1 — — A(fco-R) cos 29 


( 8 ) 


while the eigenenergy is E/{hft) = —1/2 — A^/8. 

Expanding the eigenfunctions in the basis of the plane- 
wave states, we have also solved the eigenvalue prob¬ 
lem numerically. The eigenfunctions of lowest energy 
have the general form 'b-|- = ‘- 2 m cos 2ni9 and 'bj. = 

J2m d 2 m+i cos(2m -I- 1)9, or vice versa, i.e., with Ti ex¬ 
changed with From these expressions some generic 
features follow for dt-j- and dtj,, which are also valid for the 
solutions found in Eqs. © and ®: (i) These solutions 
are even (the Hamiltonian is invariant under the trans¬ 
formation 9^—9 and thus the eigenstates are parity 
eigenstates), (ii) The density nj. has nodes at 0 = 7r/2 
and 37r/2, while at the same points ni has maxima, or 
vice versa, (hi) The periodicity of d/.|. is tt and of 
it is 27r, or vice versa, (iv) Finally, the density of both 
components has a periodicity of tt. 

Turning to the degeneracy of these solutions, the 
lowest-energy eigenstates are even. For all current ex¬ 
periments, the laser beam has a wavelength of order 1 
/im, which means that fco ^ ^tt (/im)“^. For a typical 
value of i? ^ 10 /im, kgH is thus at least 10, or larger. 
In the typical limiting case fcgi? 1, the system is in 
the striped-like phase, the density variations are located 
around opposite ends of the ring, while the density effec¬ 
tively vanishes elsewhere [see, e.g., Fig. 1(e)]. There 
is thus a corresponding antisymmetric solution which 
becomes nearly degenerate with the symmetric one in 
this limit (with a corresponding two-fold degeneracy). 
In addition, the transformation dt-j- —>• v[/j^ ——\If.|- 

and H —>■ — n leaves the eigenvalue problem unaffected. 
This gives rise to another two-fold (near) degeneracy for 
“small” H, which is exact for H = 0. 

The physical picture that emerges from the above cal¬ 
culation is that for kgR 1 (and any value of fi), i.e., 
for sufficiently small values of kg and/or of R, as well 
as for hko/{MilR) <C 1, i.e., for sufficiently small kg 
and/or large values of and/or large values of R, the 
density variation of both components is sinusoidal. As 
kgR increases (with all the other parameters fixed), the 
two components start to localize around 9 — 'k/2 and 
9 = 37r/2, with a pronounced density variation around 
these two points, which is the analogue of the striped 
phase. This transition is seen in Fig. 1, where we plot the 
one-dimensional density of the two components for an in¬ 
creasing value of kg, kgR — 1.0 (a), 1.5 (b) and 1.6 (c), 
2.8 (d), and 6 (e), and for a fixed value of h/{Ml}R^) = 1. 
The population imbalance also decreases with increasing 
kg, since the term in the Hamiltonian that is proportional 
to n becomes smaller compared to the one that is pro¬ 
portional to kg. It is also remarkable that an increase of 
the value of kgR of less than a factor of two results in a 







FIG. 1: One-dimensional density of the two pseudospin com¬ 
ponents, “up” (dotted, lighter curve) and “down” (solid, 
darker curve), within the quasi-one-dimensional model, for 
fixed values of R and fl, with fi/{MQ,R?) = 1, and for an 
increasing value of kg, kgR = 1.0 (a), 1.5 (b), 1.6 (c), 2.8 (d), 
and 6 (e). 
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rather dramatic change in the density distribution of the 
two components, as seen, e.g., between Figs. 1 (b) and 
1 (c). As we show also below, this indicates a very rich 
structure, which is a generic feature of this problem. 

In the first panel (a) of Fig. 2 we have plotted the 
density exp[-(p - i?)Vao]/(V^aoi?), 

where are the eigensolutions of Eqs. ([6]). The pa¬ 
rameters we have chosen are the ones which correspond 
to the results shown in the lower panels of Fig. 2, namely 
h/{MVtR^) = 1/16, n/w = 1, and k^R — 0.4 (first col¬ 
umn), Ti/{MQ.R?) = 1/16, n/w = 1, and k^R — 10 (sec¬ 
ond column), and finally fi/{M^R?) = 1/160, fl/w = 10, 
and k^R = 10 (third column). Panel (a) of Fig. 2 also 
serves as a reference to the lower panels which show so¬ 
lutions of the two-dimensional problem, as we discuss in 
the following section. 


IV. EIGENVALUE PROBLEM IN THE CASE OF 
AN ANNULAR POTENTIAL 

We now investigate the effect of the finiteness of the 
width of the annulus, still in the absence of interac¬ 
tions. In this case, the width of the annulus (the os¬ 
cillator length oo) provides an extra length scale, which 
is a fraction of R. For a narrow annulus (i? ^ oq) 
one goes back to the problem of quasi-one-dimensional 
motion discussed above. When oq ~ i? then the prob¬ 
lem reduces to that of a harmonic trapping potential 
[H, m, in, n HO) El) EO-fllj with a small “hole” in the 
center of the trap. The most interesting case is thus the 
intermediate one when ao ^ R- 

We show in panel (b) of Fig. 2 the two-dimensional 
density distribution of the two components that we have 
obtained numerically using the method of imaginary-time 
propagation (for the details of this calculation see the Ap¬ 
pendix), considering an annular potential with R/a^ = 4 
and some representative values of k^ and Q. In the first 
column, of “small” ko, koR = 0.4, there is a large pop¬ 
ulation imbalance between the two components - here 
h^l/Efi = 200. While the “up” component has a pro¬ 
nounced axial asymmetry, the dominant “down” compo¬ 
nent also lacks axial symmetry with the maximum of the 
density being at 9 — 0 and tt (which is very weak, and 
is hardly visible in the plot). In the second column with 
a “large” value of ko, koR = 10, both components show 
a striped phase and have an almost equal population; 
here Tifl/En = 8/25. Finally the phase shown in the 
third column with il/uj = 10 and koR = 10 resembles 
in a sense the phase with a sinusoidal density distribu¬ 
tion found within the quasi-one-dimensional model, with 
a large population imbalance (here hU/Eji = 16/5). 


V. EFFECT OF THE INTERACTIONS 

Having established a rather complete picture about 
the single-particle eigenvalue problem, we now turn to 



FIG. 2: Two-dimensional density of the two pseudospin com¬ 
ponents, “up” (higher) and “down” (lower). Panel (a) shows 
the result of the quasi-one-dimensional model with a Gaussian 
transverse profile. Panel (b) shows the density in an annular 
potential and in the absence of interactions. Panel (c) shows 
the result in an annular potential in the presence of interac¬ 
tions for g < g-fi and panel (d) for g > g-[i. The color scale is 
not the same in all the plots, since in some of them there is a 
large population imbalance. 


the effect of the interactions. In a homogeneous system 
when g-fi > g the striped phase (where also the den¬ 
sity minima/maxima of the one component coincide with 
the maxima/minima of the other) is energetically more 
favourable, since this configuration minimizes the over¬ 
lap between the two components [due to the last term 
in Eq. © below]. This may be seen if one writes the 
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FIG. 3: Two-dimensional density (black contour lines) and 
phase of the two pseudospin components, “up” (left) and 
“down” (right) for the data of the third column of Fig. 2 (d). 
The color scale shows the phase with blue as zero and red as 
27r. 


interaction energy in the following form 

-Bint = ( 5 / 2 ) J n^ir) dr + - g) j nt(r)n^(r) dr, (9) 

where n(r) = n-|-(r)-|-n 4 ,(r) is the total density. According 
to the numerical results that we have obtained using the 
method of imaginary-time propagation, when g^^^ > g the 
interactions “combine” with the (spatially dependent) 
term of the spin-orbit Hamiltonian that is proportional 
to fco (and favours the formation of stripes, as we argued 
earlier). Both terms make it energetically favourable for 
the system to form a striped-like phase (similar to the 
one that we described in the case of an annular poten¬ 
tial with spin-orbit coupling, but without interactions). 
Panel (c) of Fig. 2 shows n-j- and nj,, for = 1/2, and 
Ng/[27:Ra^huj) = 25/(47r), or Ng/{271 Ra^Ea) = 1250/7r 
in the first column and Ng/{27rRaoEii) = 2 / 7 r in the 
second and in the third column. While in the first two 
columns the density is slightly affected by the interac¬ 
tions [i.e., as compared to the density in Fig. 2 (b)], in 
the third column the interactions expand the two com¬ 
ponents more around the annulus. 

In panel (d) of Fig. 2 we have chosen gjg^i = 2. The 
first column indicates that the density is not affected 
drastically by the interactions. On the other hand, in 
the second column we see that the density is essentially 
axially symmetric, in contrast to the striped-like phase, 
which is seen in the absence of interactions in the second 
column of panel (b). Thus, this phase is strongly affected 
by the interactions. 

At this point it is interesting to make contact between 
the results of Fig. 2 and those of no trapping potential. 
As we mentioned also above, in the absence of interac¬ 
tions when the ratio hfl/{4:Efi) is smaller or larger than 
unity, the system is in the striped, or in the homogeneous 
phase, respectively [l^. The actual value of this ratio is 
50 (first column), 0.08 (second column), and 0.8 (third 
column). According to the above criterion, the first and 
the second columns of Fig. 2 (b) are in the homogeneous 
and in the striped phases, respectively, which is consis¬ 
tent with our data. For the interacting problem, the 
critical value for the ratio hfl/{4Efi) is 0.19/4 = 0.0475 
in this case [l^, [s^, Therefore, the first and the 


third columns of Fig. 2(c) and (d) are clearly in the homo¬ 
geneous phase, which again is consistent with our data. 

In the third column of Fig. 2 (d) we show a solution 
of nonzero circulation, which is only possible in the pres¬ 
ence of interactions and for a wide enough annulus. For 
the same data, Fig. 3 shows the phase of the two order 
parameters along with the corresponding density. As 
seen in this plot, there are two phase singularities in 
the “down” component and three in the “up”. In the 
“down” component the corresponding vortices reside in 
the central region of exponentially small density, while 
in the “up” component one of the three is in the cen¬ 
tral region and two in the region of nonzero density (for 
this reason the deviation of the density of the “up” com¬ 
ponent from axial symmetry is more pronounced than 
that of the “down” component). For the same reason, 
while all the other plots of Fig. 2 have a mirror symme¬ 
try with respect to the x axis, this symmetry is broken. 
We should also mention that the sense of circulation of 
the state of the third column of Fig. 2 (d) is determined 
by the initial condition that is given in the (iterative) 
method of imaginary-time propagation. Furthermore, 
with a “mean” density no = N/{277Ra^), for the heal¬ 
ing length ^0 that corre spon ds to t he coupling constant 
<?! ^o/oq = \/?iw/ {2nog) = y^27r/25 ~ 1/2. Therefore, ^0 
is roughly one half of the width of the annulus, since Og 
sets this scale. Clearly these length scales play a crucial 
role, since the vortex states we have found have to “fit” 
within the annulus and this becomes possible provided 
that the width of the annulus is wide enough, i.e., it is 
comparable, or larger than the coherence length ^o- 

The vortex states we have found in Fig. 3 involve two 
large parameters: Rabi coupling il/w = 10 and Raman 
wavenumber koR = 10 , as well and the important effect 
of interactions with g > g-^^. Each of these plays an 
important role, since the vortex states appear only in 
the right column of Fig. 2(d). As seen in Fig. 3, there are 
persistent currents in both the upper component (arising 
from the phase singularity in the central region) and the 
lower component (arising from the two phase singularities 
in the central region). In addition, the upper component 
has two vortices in the annular region with finite number 
density. Such states [H, [dO, El, El| have the potential to 
serve as cold-atom analogs of ring currents and definitely 
merit further more detailed investigation. 


VI. CONCLUSIONS 

Spin-orbit coupled Bose-Einstein-condensed atoms 
confined in an annular potential show a variety of phases 
as a result of the (single-particle effect of the) spin-orbit 
coupling and of the (many-body effect of the) coupling 
between the two atoms, while the non-trivial topology 
of the annular potential is another novel aspect of the 
problem that we have solved. 

Figure 1 summarizes the results for the solutions of 
the eigenvalue problem in the case of a very narrow an- 
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nulus, as fco increases (from top to bottom). The homo¬ 
geneous phase that is present in an untrapped system is 
replaced by a sinusoidal density variation in the case of 
a very narrow annulus, which evolves continuously into 
a striped-like phase as fco increases. As the width of the 
annulus increases, this picture persists qualitatively, as 
seen in Fig. 2, where interactions have also been included. 
Depending on the relative strength of g and the in¬ 
teractions either favor the striped phase, or suppress it, 
while they may also give rise to nonlinear solutions with 
a nonzero circulation that require nonzero interactions. 
The very interesting implication of these solutions is that 
in a spin-orbit coupled system vortex states might spon¬ 
taneously be created in one or both components, some¬ 
times in the physical region with nonzero density and 
sometimes in the central region with zero density (these 
latter vortices lead to quantized circulation around the 
annulus). 

While our study does not in any way exhaust all the 
possible phases, it gives a sense of the richness of this 
problem. A natural future project is a more system¬ 
atic study of the solutions presented in Fig. 2 (d) and 
Fig. 3. Finally, the dimensional reduction we have per¬ 
formed provides a general method, which (with the trivial 
inclusion of interactions) may be useful in investigating 
questions related with e.g., nonlinear, solitary-wave so¬ 
lutions in quasi-one-dimensional spin-orbit coupled sys¬ 
tems. 
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Appendix: Numerical method 

The numerical method used to solve the two coupled 
Gross-Pitaevskii-like equations for the two pseudospin 
components of the order parameter, given by Eq. o 
of the main text, is based on a fourth-order split-step 
Fourier method within an imaginary-time propagation 
technique. The details of the method were previously 
discussed in Ref. . Therefore here we give only a brief 
overview of the method. 

The starting point, which forms the basis of the 
imaginary-time propagation method, is to consider the 
time-dependent Schrodinger equation in imaginary time, 
i.e., T = —it (in the expressions below we use atomic 
units for convenience) 

^ (A.l) 

where H is the Hamiltonian, and to propagate the wave 


function ip under the action of a time evolution operator 
exp(-TiJ), 

^(r) oc exp{—TH)'tjj{0). (A.2) 

If the initial state ip{0) is expanded in the eigenstates (pi 
of the Hamiltonian, with eigenvalues £i , then 

V'(r) oc ^ Ci(pi exp(-rej). (A.3) 

i 

Therefore, the time-evolution operator leads to an expo¬ 
nential decay of the eigenstates, where the one with the 
lowest energy eo decays with the slowest rate and ipir) 
eventually converges to the ground state (pa as r —>■ oo. 

To evaluate the time evolution given in Eq. (IA.3I) . we 
use a forward fourth-order algorithm called split-step 
Eourier method: 

•!/'(Ar) = x 

X g-(l/2)ArTg-(l/6)Ary(O)^(0)^ (^.4) 

In our study V = H — K, where K is the kinetic energy, 
serves as an effective potential and it consists of the ex¬ 
ternal trapping potential, the interaction and the spin- 
orbit coupling terms. The mid-point effective potential 
in Eq. (jA.4|) is given by 

V = V + ^[V,[T,V]]. (A.5) 

The order parameter is discretized on a square mesh of 
2 ^ points {Af is an integer) in both x and y directions 
with separation Al and its Eourier transform is computed 
using the discretized fast Fourier transform (FFT). 

In our calculations we have used 256^ grid points with 
Al = 0.1, where we have observed that larger grid sizes 
do not have a significant effect on the obtained results. 
The time step At had to be chosen small enough to en¬ 
sure that the Hamiltonian of the system does not change 
dramatically between the iterations. The optimum value 
of At, which satisfies this condition and also provides a 
good convergence in our study, has been determined as 
0.01. On the other hand, we have observed that the use 
of smaller time steps does not improve the accuracy of 
our results. 

The time evolution procedure requires an initial condi¬ 
tion for the wave function that is propagated in imaginary 
time until a steady-state solution with the (local or ab¬ 
solute) minimum of energy is reached after a sufficiently 
large number of iterations. For example, to obtain the so¬ 
lution with nonzero circulation given in the third column 
of Fig. 2 (d) and in Fig. 3, we have performed calcu¬ 
lations using ten different initial states. We have found 
that all the results exhibit density distributions with a 
nonzero circulation and the result presented in Fig. 2 (d) 
and in Fig. 3 is the energetically most favorable state. 
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